The X-Philes (2nd Revision)/The X-Philes Number 1 (1995).iso/xphiles/hp48hor1/bessel.src

Text File | 1991-10-19 | 18KB | 672 lines

%%HP: T(3)A(D)F(.);
DIR
@ BESSEL by Mark A. Ordal
@ ----------------------------------------------------------------------
@
@                J (x), Y (x), I (x), and K (x)
@                 n      n      n          n
@
@ ----------------------------------------------------------------------
@  Integer order Bessel Functions of the first and second kinds and
@  integer order Modified Bessel Functions of the first and second kinds
@  for non-negative integer order and non-negative real argument.
@ ----------------------------------------------------------------------
@ HP48 programs by Dr. Mark A. Ordal, North Dakota State University,
@ Physics Department, Fargo, ND, 58105
@ NU123952@NDSUVM1
@ revision of April 25, 1991
@ ----------------------------------------------------------------------
@ This is the ASCII listing (translate 3 mode) of a directory that
@ contains the following:
@
@    Programs:
@        Jnx, NRJN, ASJ1, and ASJ0      (Bessel: first kind)
@        Ynx, NRYN, ASY1, and ASY0      (Bessel: second kind)
@        Inx, NRIN, ASI1, and ASI0      (Modified Bessel: first kind)
@        Knx, NRKN, ASK1, and ASK0      (Modified Bessel: second kind)
@
@    Subroutines: JYIK, JY1, and JY0
@
@ NOTE: These programs do not prevent you from entering an invalid order
@       or argument.
@
@ When named BESSINT, running the Bytes command on the name of this
@ directory returns:  Checksum  #18596 and 4558 bytes.
@ ----------------------------------------------------------------------
@ Changes from previous version:
@   1) Programs ASJ0, ASJ1, ASY0, and ASY1 preserve the state of system
@      flags (e.g., the DEG or RAD modes).
@
@   2) Subroutine JYIK was made 4.5 bytes smaller by using a START-NEXT
@      LOOP instead of a FOR-NEXT loop.
@ ----------------------------------------------------------------------
Jnx
@ ----------------------------------------------------------------------
@ Bessel Functions of the FIRST kind and integer order
@
@ To use:
@ Level 2 of stack: order (a nonnegative integer)
@ Level 1 of stack: argument (any nonnegative real number)
@
@ When named Jnx, the BYTES command returns #51151 and 123.5 bytes
@ ----------------------------------------------------------------------
  Jnx
    \<< \-> n x
      \<<
        CASE n 0
SAME
          THEN x
ASJ0
          END n 1
SAME
          THEN x
ASJ1
          END n x
NRJN
        END
      \>>
    \>>
  Ynx
@ ----------------------------------------------------------------------
@ Bessel Functions of the SECOND kind and integer order
@
@ To use:
@ Level 2 of stack: order (a nonnegative integer)
@ Level 1 of stack: argument (any positive real number)
@
@ Remember that Yn(x) is infinite for x=0
@
@ When named Ynx, the BYTES command returns #18975 and 123.5 bytes
@ ----------------------------------------------------------------------
    \<< \-> n x
      \<<
        CASE n 0
SAME
          THEN x
ASY0
          END n 1
SAME
          THEN x
ASY1
          END n x
NRYN
        END
      \>>
    \>>
  Inx
@ ----------------------------------------------------------------------
@ Modified Bessel Functions of the FIRST kind and integer order
@
@ To use:
@ Level 2 of stack: order (a nonnegative integer)
@ Level 1 of stack: argument (any nonnegative real number)
@
@ When named Inx, the BYTES command returns #30214d and 123.5 bytes
@ ----------------------------------------------------------------------
    \<< \-> n x
      \<<
        CASE n 0
SAME
          THEN x
ASI0
          END n 1
SAME
          THEN x
ASI1
          END n x
NRIN
        END
      \>>
    \>>
  Knx
@ ----------------------------------------------------------------------
@ Modifed Bessel Functions of the SECOND kind and integer order
@
@ To use:
@ Level 2 of stack: order (a nonnegative integer)
@ Level 1 of stack: argument (any positive real number)
@
@ Remember that Kn(x) is infinite for x=0
@
@ When named Knx, the BYTES command returns #20615d and 123.5 bytes
@ ----------------------------------------------------------------------
    \<< \-> n x
      \<<
        CASE n 0
SAME
          THEN x
ASK0
          END n 1
SAME
          THEN x
ASK1
          END n x
NRKN
        END
      \>>
    \>>
  NRJN
@ ----------------------------------------------------------------------
@ Bessel Functions of the FIRST kind and integer order n > 1
@ calculated using the RPL equivalent of the "Numerical Recipes"
@ program BESSJ.  Needs the equivalent of the "Numerical Recipes"
@ programs BESSJ0 and BESSJ1.  In this case, programs ASJ0 and ASJ1
@ are used instead -- based on Abramowitz and Stegun Eqs 9.4.1 and
@ 9.4.3 for order zero and Eqs 9.4.4 and 9.4.6 for order one.
@
@ The programs ASJ0 and ASJ1 based on the Abramowitz and Stegun
@ equations are shorter and faster than RPL equivalents of BESSJ0
@ and BESSJ1.)
@
@ To use:
@ Level 2 of stack: order (a positive integer > 1)
@ Level 1 of stack: argument (any nonnegative real number)
@
@ The "Numerical Recipes" program variables TOX, SUM, BIGNO, BIGNI,
@ and BESSJ have been replaced by t, s, u, d, and sj, respectively.
@ Variable IACC has been replaced by its value of 40.
@
@ When named NRJN, the BYTES command returns #18067d and 657 bytes
@
@ To adapt for use on the HP28S, replace the commands 'm' DECR with
@ the equivalent sequence m 1 - 'm' STO m
@ ----------------------------------------------------------------------
    \<< 10000000000
.0000000001 0 0 0 0
\-> n x u d t m s sj
      \<< 2 x / 't'
STO
        IF x n >
        THEN x ASJ0
x ASJ1 1 n 1 -
          FOR j j t
* OVER * 3 PICK -
ROT DROP
          NEXT SWAP
DROP
        ELSE 40 n *
\v/ IP n + 2 / IP 2 *
'm' STO 0 1 m
          DO t *
OVER * 3 PICK - ROT
DROP
            IF DUP
ABS u >
            THEN d
* SWAP d * SWAP d
DUP 'sj' STO* 's'
STO*
            END
            IF m n
SAME
            THEN
OVER 'sj' STO
            END 'm'
DECR t * OVER * 3
PICK - ROT DROP
            IF DUP
ABS u >
            THEN d
* SWAP d * SWAP d
DUP 'sj' STO* 's'
STO*
            END
            IF m n
SAME
            THEN
OVER 'sj' STO
            END DUP
's' STO+ 'm' DECR
          UNTIL m 1
<
          END DROP
SWAP DROP NEG s 2 *
+ sj SWAP /
        END
      \>>
    \>>
  NRYN
@ ----------------------------------------------------------------------
@ Bessel Functions of the SECOND kind and integer order n > 1
@ calculated using the RPL equivalent of the "Numerical Recipes"
@ program BESSY.  Needs the equivalent of the "Numerical Recipes"
@ programs BESSY0 and BESSY1.  In this case, programs ASY0 and ASY1
@ are used instead -- based on Abramowitz and Stegun Eqs 9.4.1 and
@ 9.4.3 for order zero and Eqs 9.4.4 and 9.4.6 for order one.
@
@ The programs ASY0 and ASY1 based on the Abramowitz and Stegun
@ equations are shorter and faster than RPL equivalents of BESSY0
@ and BESSY1.)
@
@ To use:
@ Level 2 of stack: order (a positive integer > 1)
@ Level 1 of stack: argument (any positive real number)
@
@ The "Numerical Recipes" program variable TOX has been replaced by t.
@ Variable IACC has been replaced by its value of 40.
@
@ When named NRYN, the BYTES command returns #16760 and 143 bytes
@ ----------------------------------------------------------------------
    \<< 0 \-> n x t
      \<< 2 x / 't'
STO x ASY0 x ASY1 1
n 1 -
        FOR j j t *
OVER * 3 PICK - ROT
DROP
        NEXT SWAP
DROP
      \>>
    \>>
  NRIN
@ ----------------------------------------------------------------------
@ Modified Bessel Functions of the FIRST kind and integer order n > 1
@ calculated using the RPL equivalent of the "Numerical Recipes"
@ program BESSI.  Needs the equivalent of the "Numerical Recipes"
@ programs BESSI0 and BESSI1.  These programs (renamed ASJ0 and ASJ1)
@ are based on Abramowitz and Stegun Eqs 9.8.1 and 9.8.2 for order zero
@ and Eqs 9.8.3 and 9.8.4 for order one.
@
@ To use:
@ Level 2 of stack: order (a positive integer > 1)
@ Level 1 of stack: argument (any nonnegative real number)
@
@ The "Numerical Recipes" program variables TOX, BIGNO, BIGNI,
@ and BESSI have been replaced by t, u, d, and si, respectively.
@ Variable IACC has been replaced by its value of 40.
@
@ When named NRIN, the BYTES command returns #9623d and 315 bytes
@ ----------------------------------------------------------------------
    \<< 10000000000
.0000000001 0 0 \-> n
x u d t si
      \<< 2 x / 't'
STO 0 1 40 n * \v/ IP
n + 2 * 1
        FOR j j t *
OVER * 3 PICK + ROT
DROP
          IF DUP
ABS u >
          THEN d *
SWAP d * SWAP d
'si' STO*
          END
          IF j n
SAME
          THEN OVER
'si' STO
          END -1
        STEP SWAP
DROP si SWAP / x
ASI0 *
      \>>
    \>>
  NRKN
@ ----------------------------------------------------------------------
@ Modifed Bessel Functions of the SECOND kind and integer order n > 1
@ calculated using the RPL equivalent of the "Numerical Recipes"
@ program BESSK.  Needs the equivalent of the "Numerical Recipes"
@ programs BESSK0 and BESSK1.  These programs (renamed ASK0 and ASK1)
@ are based on Abramowitz and Stegun Eqs 9.8.5 and 9.8.6 for order zero
@ and Eqs 9.8.7 and 9.8.8 for order one.
@
@ To use:
@ Level 2 of stack: order (a positive integer > 1)
@ Level 1 of stack: argument (any positive real number)
@
@ The "Numerical Recipes" program variables TOX, BIGNO, BIGNI,
@ and BESSK have been replaced by t, u, d, and sk, respectively.
@
@ When named NRKN, the BYTES command returns #20286d and 181 bytes
@ ----------------------------------------------------------------------
    \<< 10000000000
.0000000001 0 0 \-> n
x u d t sk
      \<< 2 x / 't'
STO x ASK0 x ASK1 1
n 1 -
        FOR j j t *
OVER * 3 PICK + ROT
DROP
        NEXT SWAP
DROP
      \>>
    \>>
  ASJ1
@ ----------------------------------------------------------------------
@ Bessel Functions of t4e FIRST kind and order one calculated
@ using Abramowitz and Stegun Eqs 9.4.4 and 9.4.6
@
@ This program is slightly shorter and faster than the RPL equivalent
@ of the "Numerical Recipes" program BESSJ1.
@
@ This program program calls programs JYIK and JY1.
@
@ To use:
@ Level 1 of stack: argument (any nonnegative real number)
@
@ When named ASJ1, the BYTES command returns #49835d and 213.5 bytes
@ ----------------------------------------------------------------------
    \<< 0 RCLF \-> x a
ff
      \<<
        IF x 3 <
        THEN .5
-.56249985
.21093573
-.03954289
.00443319
-.00031761
.00001109 x 3 / SQ
6 JYIK x *
        ELSE x JY1
RAD COS * x \v/ /
        END ff STOF
      \>>
    \>>
  ASY1
@ ----------------------------------------------------------------------
@ Bessel Functions of the SECOND kind and order one calculated
@ using Abramowitz and Stegun Eqs 9.4.4 and 9.4.6
@
@ This program is slightly shorter and faster than the RPL equivalent
@ of the "Numerical Recipes" program BESSY1.
@
@ This program program calls programs JYIK, JY1, and ASJ1.
@
@ To use:
@ Level 1 of stack: argument (any positive real number)
@
@ When named ASY1, the BYTES command returns #1027d and 270 bytes
@ ----------------------------------------------------------------------
    \<< 0 RCLF \-> x a
ff
      \<<
        IF x 3 <
        THEN
-.6366198 .2212091
2.1682709
-1.3164827 .3123951
-.0400976 .0027873
x 3 / SQ 6 JYIK x
ASJ1 x .5 * LN * x
* 2 * \pi \->NUM / + x
/
        ELSE x JY1
RAD SIN * x \v/ /
        END ff STOF
      \>>
    \>>
  ASI1
@ ----------------------------------------------------------------------
@ Modifed Bessel Functions of the FIRST kind and order one calculated
@ using Abramowitz and Stegun Eqs 9.8.3 and 9.8.4
@
@ To use:
@ Level 1 of stack: argument (any nonnegative real number)
@
@ This program program calls program JYIK
@
@ When named ASI1, the BYTES command returns #12477d and 334.5 bytes
@ ----------------------------------------------------------------------
    \<< 0 \-> x a
      \<<
        IF x ABS
3.75 <
        THEN .5
.87890594 .51498869
.15084934 .02658733
.00301532 .00032411
x 3.75 / SQ 6 JYIK
x *
        ELSE
.39894228
-.03988024
-.00362018
.00163801
-.01031555
.02282967
-.02895312
.01787654
-.00420059 3.75 x
ABS / 8 JYIK x ABS
DUP EXP SWAP \v/ / *
        END
      \>>
    \>>
  ASK1
@ ----------------------------------------------------------------------
@ Modifed Bessel Functions of the SECOND kind and order one calculated
@ using Abramowitz and Stegun Eqs 9.8.7 and 9.8.8
@
@ This program program calls programs JYIK and ASI1.
@
@ To use:
@ Level 1 of stack: argument (any positive real number)
@
@ When named ASK1, the BYTES command returns #54561d and 308 bytes
@ ----------------------------------------------------------------------
    \<< 0 \-> x a
      \<<
        IF x ABS 2
<
        THEN 1
.15443144
-.67278579
-.18156897
-.01919402
-.00110404
-.00004686 x 2 / SQ
6 JYIK x / x ASI1 x
2 / LN * +
        ELSE
1.25331414
.23498619 -.0365562
.01504268
-.00780353
.00325614
-.00068245 2 x / 6
JYIK x DUP NEG EXP
SWAP \v/ / *
        END
      \>>
    \>>
  ASJ0
@ ----------------------------------------------------------------------
@ Bessel Functions of the FIRST kind and order zero calculated
@ using Abramowitz and Stegun Eqs 9.4.1 and 9.4.3
@
@ This program is slightly shorter and faster than the RPL equivalent
@ of the "Numerical Recipes" program BESSJ0.
@ This program program calls programs JYIK and JY0.
@
@ To use:
@ Level 1 of stack: argument (any nonnegative real number)
@
@ When named ASJ0, the BYTES command returns #1504d and 198.5 bytes
@ ----------------------------------------------------------------------
    \<< 0 RCLF \-> x a
ff
      \<<
        IF x 3 <
        THEN 1
-2.2499997
1.2656208 -.3163866
.0444479 -.0039444
.00021 x 3 / SQ 6
JYIK
        ELSE x JY0
RAD COS * x \v/ /
        END ff STOF
      \>>
    \>>
  ASY0
@ ----------------------------------------------------------------------
@ Bessel Functions of the SECOND kind and order zero calculated
@ using Abramowitz and Stegun Eqs 9.4.1 and 9.4.3
@
@ This program is slightly shorter and faster than the RPL equivalent
@ of the "Numerical Recipes" program BESSY0.
@ This program program calls programs JYIK, JY0, and ASJ0.
@
@ To use:
@ Level 1 of stack: argument (any nonnegative real number)
@
@ When named ASY0, the BYTES command returns #37606d and 256 bytes
@ ----------------------------------------------------------------------
    \<< 0 RCLF \-> x a
ff
      \<<
        IF x 3 <
        THEN
.36746691 .60559366
-.74350384
.25300117
-.04261214
.00427916
-.00024846 x 3 / SQ
6 JYIK x ASJ0 x .5
* LN * 2 * \pi \->NUM /
+
        ELSE x JY0
RAD SIN * x \v/ /
        END ff STOF
      \>>
    \>>
  ASI0
@ ----------------------------------------------------------------------
@ Modifed Bessel Functions of the FIRST kind and order zero calculated
@ using Abramowitz and Stegun Eqs 9.8.1 and 9.8.2
@
@ To use:
@ Level 1 of stack: argument (any nonnegative real number)
@
@ This program program calls program JYIK
@
@ When named ASI0, the BYTES command returns #63274 and 319.5 bytes
@ ----------------------------------------------------------------------
    \<< 0 \-> x a
      \<<
        IF x ABS
3.75 <
        THEN 1
3.5156229 3.0899424
1.2067492 .2659732
.0360768 .0045813 x
3.75 / SQ 6 JYIK
        ELSE
.39894228 .01328592
.00225319
-.00157565
.00916281
-.02057706
.02635537
-.01647633
.00392377 3.75 x
ABS / 8 JYIK x ABS
DUP EXP SWAP \v/ / *
        END
      \>>
    \>>
  ASK0
@ ----------------------------------------------------------------------
@ Modifed Bessel Functions of the SECOND kind and order zero calculated
@ using Abramowitz and Stegun Eqs 9.8.5 and 9.8.6
@
@ This program program calls programs JYIK and ASI0.
@
@ To use:
@ Level 1 of stack: argument (any positive real number)
@
@ When named ASK0, the BYTES command returns #64307 and 311.5 bytes
@ ----------------------------------------------------------------------
    \<< 0 \-> x a
      \<<
        IF x ABS 2
<
        THEN
-.57721566 .4227842
.23069756 .0348859
.00262698 .0001075
.0000074 x 2 / SQ 6
JYIK x ASI0 x 2 /
LN NEG * +
        ELSE
1.25331414
-.07832358
.02189568
-.01062446
.00587872 -.0025154
.00053208 2 x / 6
JYIK x DUP NEG EXP
SWAP \v/ / *
        END
      \>>
    \>>
  JYIK
@ ----------------------------------------------------------------------
@ Subprogram for use by ASJ0, ASJ1, ASY0, ASY1, ASI0, ASI1, ASK0, and
@ ASK1
@
@ When named JYIK, the BYTES command returns #32125d and 56.5 bytes
@ ----------------------------------------------------------------------
    \<< \-> t j
      \<< 1 j
        START t * +
        NEXT
      \>>
    \>>
  JY1
@ ----------------------------------------------------------------------
@ Subprogram for use by ASJ1 and ASY1
@
@ This program program calls program JYIK
@
@ When named JY1, the BYTES command returns #54172d and 241 bytes
@ ----------------------------------------------------------------------
    \<< 0 \-> x a
      \<< 3 x / 'a'
STO .79788456
.00000156 .01659667
.00017105
-.00249511
.00113653
-.00020033 a 6 JYIK
-2.35619449
.12499612 .0000565
-.00637879
.00074348 .00079824
-.00029166 a 6 JYIK
x +
      \>>
    \>>
  JY0
@ ----------------------------------------------------------------------
@ Subprogram for use by ASJ0 and ASY0
@
@ This program program calls program JYIK
@
@ When named JY0.SUB, the BYTES command returns #15249d and 241 bytes
@ ----------------------------------------------------------------------
    \<< 0 \-> x a
      \<< 3 x / 'a'
STO .79788456
-.00000077
-.0055274
-.00009512
.00137237
-.00072805
.00014476 a 6 JYIK
-.78539816
-.04166397
-.00003954
.00262573
-.00054125
-.00029333
.00013558 a 6 JYIK
x +
      \>>
    \>>
END